Numerical simulation apparatus for time dependent schrödinger equation

ABSTRACT

To provide a numerical simulation apparatus capable of executing a numerical simulation with high speed and precision by reducing computational complexity. A numerical simulation apparatus that executes a numerical simulation using a wave function which is a solution of a time dependent Schrödinger equation includes: a real time evolution calculation unit that calculates a second wave function while evolving the second wave function from an initial time in increments of a predetermined time period, the second wave function being obtained by applying a central difference approximation in a real-space finite-difference method to a first wave function expressed using a propagator, and being expressed using a Bessel function; and a calculation result storage unit that stores a calculation result of the second wave function obtained at each time by the time evolution calculation unit while evolving the second wave function in increments of the predetermined time period.

TECHNICAL FIELD

The present invention relates to a numerical simulation apparatus thatis needed for developing electronic devices or quantum computers inconsideration of nanoscale properties, and in particular relates to anumerical simulation apparatus capable of executing a numericalsimulation with high speed and precision by reducing computationalcomplexity.

BACKGROUND ART

In recent years, active research has been conducted on nanoscaleelectronic devices, quantum computers, and the like. In the case wherean electronic device is a structure of a scale equal to or smaller thanabout 10 nm which is a mean free path of electrons, quantum effectsbecome noticeable in electrical and magnetic properties such aselectrical conductance. Thus, it is difficult to design and developthese nanoscale electronic devices and quantum computers according to anordinary method that is based on macroscale electrical and magneticproperties. Accordingly, a numerical simulation apparatus capable ofperforming a numerical simulation of dynamic behavior of electrons suchas scattering, transmission, reflection, interference, oscillation,attenuation, and excitation with high speed and precision is needed forthe design and development of these nanoscale electronic devices andquantum computers.

For instance, with miniaturization and high integration of a nanoscaleelectronic device, a conduction electron path becomes an atomic size(for example, see (a) in FIG. 1). As a result, electronic transportproperties specific to a nanoregion, such as ballistic conduction,appear. Therefore, in evaluations of leakage currents of silicon oxideand conduction properties of a quantum wire (for example, see (b) inFIG. 1), it is necessary to numerically simulate conduction propertiesusing the time dependent Schrödinger equation.

In a numerical simulation of conduction properties, the steady-stateSchrödinger equation for a scattering problem is typically used. Thisbeing so, various methods for numerically simulating the steady-stateSchrödinger equation for a scattering problem have been proposed. One ofthem is the Overbridging Boundary-Matching (OBM) method which employsthe real-space finite-difference method (for example, see Non-patentReference 1).

On the other hand, there is also a method for clarifying dynamicbehavior of electrons by directly solving the time dependent Schrödingerequation. In this case, the time dependent Schrödinger equation definedby the following Expression (1) is used.

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 1} \rbrack & \; \\{{{\mathbb{i}}\frac{\partial}{\partial t}{\Psi( {r,t} )}} = {H\;{\Psi( {r,t} )}}} & (1)\end{matrix}$

H is a Hamiltonian, and is given by the following Expression (2). Δ is aLaplacian, and is given by the following Expression (3). V(r, t) is apotential of an electron. Atomic units are used as a system of units. Anelectron mass is m=1, an elementary charge is e=1, and h/2n=1, where his a Planck's constant.

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 2} \rbrack & \; \\{H = {{{- \frac{1}{2}}\Delta} + {V( {r,t} )}}} & (2) \\\lbrack {{Expression}\mspace{14mu} 3} \rbrack & \; \\{\Delta = {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} + \frac{\partial^{2}}{\partial z^{2}}}} & (3)\end{matrix}$

Furthermore, expressing a solution of the above Expression (1) by a pathintegral using the propagator K (Feynman Kernel) leads to the followingExpression (4), where Ψ(r, t₀) is a wave function at initial time t₀.

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 4} \rbrack & \; \\{{\Psi( {r,t} )} = {\int_{- \infty}^{\infty}{{K( {r,{t;r^{\prime}},t_{0}} )}{\Psi( {r^{\prime},t_{0}} )}\ {\mathbb{d}^{3}r^{\prime}}}}} & (4)\end{matrix}$

It is known that the above Expression (4) can be expressed by a pathintegral. For instance, the wave function after short time period Δt isdefined by the following Expression (5) (for example, see Non-patentReference 2).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 5} \rbrack & \; \\{{\Psi( {r,{t + {\Delta\; t}}} )} = {\sqrt{\frac{1}{2\pi\;{\mathbb{i}}\;\Delta\; t}}{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{\mathbb{i}}\frac{{({r - r^{\prime}})}^{2}}{2\Delta\; t}}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({\frac{r + r^{\prime}}{2},t})}}}\ {\Psi( {r^{\prime},t} )}{\mathbb{d}^{3}r^{\prime}}}}}} & (5)\end{matrix}$

-   Non-patent Reference 1: K. Hirose, T. Ono et al., First-Principles    Calculations in Real-Space Formalism, Imperial College Press, London    (2005)-   Non-patent Reference 2: R. P. Feynman and A. R. Hibbs, Quantum    Mechanics and Path Integrals, translated by Kazuo Kitahara,    McGraw-Hill (July 1990)

SUMMARY OF THE INVENTION Problems that Invention is to Solve

With the conventional calculation methods, however, only small modelcalculations and qualitative evaluations have been attempted so far, dueto their heavy computational burdens. For example, in a numericalsimulation using the above Expression (5), enormous computationalcomplexity is required for a process of sequentially executingintegration and finding a correct solution after time t. Therefore, onlysmall models can be subject to the numerical simulation using the aboveExpression (5). Thus, the conventional techniques have failed to providea sufficient design tool for concrete electronic devices or quantumcomputers of a three-dimensional structure composed of atoms.

The present invention has been developed in view of the above problem,and has an object of providing a numerical simulation apparatus capableof executing a numerical simulation with high speed and precision byreducing computational complexity.

Means to Solve the Problems

To achieve the stated object, the numerical simulation apparatusaccording to the present invention is (a) a numerical simulationapparatus that executes a numerical simulation using a wave functionwhich is a solution of a time dependent Schrödinger equation, thenumerical simulation apparatus including: (a1) a time evolutioncalculation unit that calculates a second wave function while evolvingthe second wave function from an initial time in increments of apredetermined time period, the second wave function being obtained byapplying a central difference approximation in a real-spacefinite-difference method to a first wave function expressed using apropagator, and being expressed using a Bessel function; and (a2) acalculation result storage unit that stores a calculation result of thesecond wave function obtained at each time by the time evolutioncalculation unit while evolving the second wave function in incrementsof the predetermined time period.

According to this, numerical calculations of time evolution of a wavefunction according to the time dependent Schrödinger equation can beperformed with high speed and precision.

Moreover, (b) the numerical simulation apparatus may further include:(b1) an initial wave function storage unit that stores a pulse functionas the second wave function at the initial time; and (b2) a frequencyanalysis unit that calculates a steady-state solution of a scatteringproblem for predetermined energy, by performing a frequency analysis ofthe second wave function for the predetermined energy using thecalculation result of the second wave function obtained by the timeevolution calculation unit from the pulse function stored in the initialwave function storage unit.

According to this, a steady-state solution of a scattering problem canbe obtained with high speed and precision by numerical calculations.

Moreover, (c) the time evolution calculation unit may calculateimaginary time evolution of the second wave function, using imaginarytime that is obtained by multiplying time by an imaginary number.

According to this, numerical calculations of imaginary time evolution ofa wave function used when determining a ground state and an excitedstate can be performed with high speed and precision.

Moreover, (d) the numerical simulation apparatus may further include anabsorbing boundary condition storage unit that stores an absorptioncoefficient weighted so as to absorb an influence that is exerted on amodel by the wave function in a boundary region of the model, the modelbeing subject to the numerical simulation, wherein the time evolutioncalculation unit calculates the time evolution of the second wavefunction using the absorption coefficient stored in the absorbingboundary condition storage unit.

According to this, a large calculation region needed when setting aninfinite boundary condition can be significantly reduced.

It should be noted that the present invention can be realized not onlyas the numerical simulation apparatus, but also as a numericalsimulation method for controlling the numerical simulation apparatus, anumerical simulation program for causing one or more computer systemsand the like to execute the numerical simulation method, a recordingmedium on which the numerical simulation program is recorded, and thelike. The present invention can also be realized as a Large ScaleIntegration (LSI) including the functions of the numerical simulationapparatus, an Intellectual Property (IP) core for implementing thesefunctions on a programmable logic device such as a Field ProgrammableGate Array (FPGA) or a Complex Programmable Logic Device (CPLD), arecording medium on which the IP core is recorded, and the like.

Effects of the Invention

According to the present invention, by employing the wave functionexpressed using the Bessel function, numerical calculations of timeevolution of a wave function according to the time dependent Schrödingerequation can be performed with high speed and precision. In addition, asteady-state solution of a scattering problem can be obtained with highspeed and precision by numerical calculations. Further, numericalcalculations of imaginary time evolution of a wave function used whendetermining a ground state and an excited state can be performed withhigh speed and precision. Moreover, a large calculation region neededwhen setting an infinite boundary condition can be significantlyreduced.

Furthermore, by applying the present invention to the density functionaltheory which is a dominant method of electronic state calculationstoday, a numerical solution of the Kohn-Sham equation can be obtainedwith high speed. This makes it possible to numerically simulate dynamicbehavior of electrons for an actual model. Thus, the present inventioncan be used for design of nanoscale electronic devices and quantumcomputers.

BRIEF DESCRIPTION OF DRAWINGS

In FIG. 1, (a) shows an overview of a nanoscale electronic device, and(b) shows an overview of a quantum wire.

FIG. 2 shows a structure of a numerical simulation apparatus in anembodiment according to the present invention.

FIG. 3 shows a numerical simulation process executed in the numericalsimulation apparatus in the embodiment according to the presentinvention.

FIG. 4 shows a value of a Bessel function and a real part of anintegrand.

FIG. 5 shows a setting example of a coefficient by which a wave functionis multiplied in an absorbing boundary condition.

FIG. 6 shows a one-dimensional scattering problem by a potentialbarrier.

FIG. 7 shows a result of analyzing a frequency of the wave functionevolved with time and a result of calculating transmittance to thepotential barrier.

FIG. 8 shows a two-dimensional wire model.

FIG. 9 shows a result of analyzing a frequency in the two-dimensionalwire model and a result of calculating conductance.

FIG. 10 shows a scattering wave function obtained by applying anabsorbing boundary condition.

FIG. 11 shows calculation cost evaluations.

NUMERICAL REFERENCES

-   -   100 Numerical simulation apparatus    -   101 Calculation condition setting unit    -   102 Calculation condition storage unit    -   103 Initial value setting unit    -   104 Initial value storage unit    -   105 Real time evolution judgment unit    -   106 Real time evolution calculation unit    -   107 Imaginary time evolution calculation unit    -   108 Calculation result storage unit    -   109 Frequency analysis unit    -   110 Analysis result storage unit    -   111 Physical property calculation unit    -   112 Calculation result storing unit    -   113 Calculation result output unit

DETAILED DESCRIPTION OF THE INVENTION Embodiment

An embodiment according to the present invention is described below,with reference to drawings.

A numerical simulation apparatus in this embodiment has the followingfeatures (a) to (d).

(a) A numerical simulation apparatus that executes a numericalsimulation using a wave function which is a solution of a time dependentSchrödinger equation, includes: (a1) a time evolution calculationfunction that calculates a second wave function while evolving thesecond wave function from an initial time in increments of apredetermined time period, the second wave function being obtained byapplying a central difference approximation in a real-spacefinite-difference method to a first wave function expressed using apropagator, and being expressed using a Bessel function; and (a2) acalculation result storage function that stores a calculation result ofthe second wave function obtained at each time by the time evolutioncalculation function while evolving the second wave function inincrements of the predetermined time period.

(b) The numerical simulation apparatus includes: (b1) an initial wavefunction storage function that stores a pulse function as the secondwave function at the initial time; and (b2) a frequency analysisfunction that calculates a steady-state solution of a scattering problemfor predetermined energy, by performing a frequency analysis of thesecond wave function for the predetermined energy using the calculationresult of the second wave function obtained by the time evolutioncalculation function from the pulse function stored in the initial wavefunction storage function.

(c) The time evolution calculation function calculates imaginary timeevolution of the second wave function, using imaginary time that isobtained by multiplying time by an imaginary number.

(d) The numerical simulation apparatus includes an absorbing boundarycondition storage function that stores an absorption coefficientweighted so as to absorb an influence that is exerted on a model by thewave function in a boundary region of the model, the model being subjectto the numerical simulation, wherein the time evolution calculationfunction calculates the time evolution of the second wave function usingthe absorption coefficient stored in the absorbing boundary conditionstorage function.

On the basis of the above, the numerical simulation apparatus in thisembodiment is described below.

FIG. 2 shows a structure of the numerical simulation apparatus in thisembodiment. As shown in FIG. 2, a numerical simulation apparatus 100includes a calculation condition setting unit 101, a calculationcondition storage unit 102, an initial value setting unit 103, aninitial value storage unit 104, a real time evolution judgment unit 105,a real time evolution calculation unit 106, an imaginary time evolutioncalculation unit 107, a calculation result storage unit 108, a frequencyanalysis unit 109, an analysis result storage unit 110, a physicalproperty calculation unit 111, a calculation result storing unit 112, acalculation result output unit 113, and the like.

The calculation condition setting unit 101 receives calculationcondition data from an operator (not illustrated) via an input device(not illustrated) or reads calculation condition data from a file (notillustrated) containing calculation condition data, when executing anumerical simulation. The calculation condition setting unit 101 outputsthe calculation condition data received from the operator or thecalculation condition data read from the file, to the calculationcondition storage unit 102. Here, the calculation condition data mayinclude a space size (L_(x), L_(y), L_(z)), a space step size (h_(x),h_(y), h_(z)), a time step size (Δt), a calculation time (T), and so on.

The calculation condition storage unit 102 stores the calculationcondition data outputted from the calculation condition setting unit101.

The initial value setting unit 103 receives an initial value from theoperator (not illustrated) via the input device (not illustrated) orreads an initial value from a file (not illustrated) containing aninitial value, when executing the numerical simulation. The initialvalue setting unit 103 outputs the initial value received from theoperator or the initial value read from the file, to the initial valuestorage unit 104. Here, the initial value may include an initial wavefunction (X(r_(j), t₀)), an initial potential (V(r_(j), t)), anabsorbing boundary condition, and so on.

The initial value storage unit 104 stores the initial value outputtedfrom the initial value setting unit 103.

The real time evolution judgment unit 105 receives instruction data fromthe operator (not illustrated) via the input device (not illustrated) orreads instruction data from a file (not illustrated) containinginstruction data, when executing the numerical simulation. The real timeevolution judgment unit 105 judges whether or not to execute a real timeevolution calculation process, from the instruction data received fromthe operator or the instruction data read from the file. When judging toexecute the real time evolution calculation process, the real timeevolution judgment unit 105 calls the real time evolution calculationunit 106 to execute the real time evolution calculation process. Whenjudging not to execute the real time evolution calculation process, thereal time evolution judgment unit 105 calls the imaginary time evolutioncalculation unit 107 to execute an imaginary time evolution calculationprocess.

The real time evolution calculation unit 106 executes the real timeevolution calculation process, when called by the real time evolutionjudgment unit 105. Here, the real time evolution calculation unit 106executes the real time evolution calculation process using thecalculation condition stored in the calculation condition storage unit102 and the initial value stored in the initial value storage unit 104.The real time evolution calculation unit 106 outputs an obtainedcalculation result to the calculation result storage unit 108.

The “real time evolution calculation process” is a process ofcalculating time evolution of a wave function, by calculating, throughthe use of the wave function at predetermined time t, the wave functionat time t+Δt which is evolved from predetermined time t by short timeperiod Δt, and further repeating this wave function calculation throughthe use of the calculated wave function.

It should be noted that a second wave function (Expression (6) shownbelow) is employed as a wave function which is a solution of the timedependent Schrödinger equation. The second wave function is obtained byapplying a central difference approximation in the real-spacefinite-difference method to a first wave function (Expression (5) shownabove) expressed by a path integral using a propagator (Feynman Kernel),and is expressed using a Bessel function.

The second wave function may also be obtained by applying the centraldifference approximation in the real-space finite-difference method to awave function expressed using a time evolution Green's function.

When the real time evolution calculation process ends, the real timeevolution calculation unit 106 further calls the frequency analysis unit109 and the physical property calculation unit 111 to execute afrequency analysis process and a physical property calculation process,respectively. When these processes end, the real time evolutioncalculation unit 106 calls the calculation result output unit 113, andends the series of processes. The imaginary time evolution calculationunit 107 executes the imaginary time evolution calculation process, whencalled by the real time evolution judgment unit 105. Here, the imaginarytime evolution calculation unit 107 executes the imaginary timeevolution calculation process using the calculation condition stored inthe calculation condition storage unit 102 and the initial value storedin the initial value storage unit 104. The imaginary time evolutioncalculation unit 107 outputs an obtained calculation result to thecalculation result storage unit 108.

The “imaginary time evolution calculation process” is a process ofcalculating imaginary time evolution of the wave function by usingimaginary time that is obtained by multiplying real time by an imaginarynumber.

When the imaginary time evolution calculation process ends, theimaginary time evolution calculation unit 107 further calls the physicalproperty calculation unit 111 to execute the physical propertycalculation process. When the physical property calculation processends, the imaginary time evolution calculation unit 107 calls thecalculation result output unit 113, and ends the series of processes.

The calculation result storage unit 108 stores the calculation resultoutputted from the real time evolution calculation unit 106. Likewise,the calculation result storage unit 108 stores the calculation resultoutputted from the imaginary time evolution calculation unit 107.

The frequency analysis unit 109 executes the frequency analysis process,when called by the real time evolution calculation unit 106. Here, thefrequency analysis unit 109 executes the frequency analysis processusing the calculation condition data stored in the calculation conditionstorage unit 102, the initial value stored in the initial value storageunit 104, and the calculation result stored in the calculation resultstorage unit 108. The frequency analysis unit 109 outputs an obtainedanalysis result to the analysis result storage unit 110.

The “frequency analysis process” is a process of calculating asteady-state solution of a scattering problem for predetermined energy,by performing a frequency analysis of the wave function for thepredetermined energy using the calculation result obtained as a resultof the real time evolution calculation process from a pulse functionthat is stored as the wave function at the initial time.

The analysis result storage unit 110 stores the analysis resultoutputted from the frequency analysis unit 109.

The physical property calculation unit 111 executes the physicalproperty calculation process, when called by the real time evolutioncalculation unit 106. Here, the physical property calculation unit 111executes the physical property calculation process using the calculationcondition data stored in the calculation condition storage unit 102, theinitial value stored in the initial value storage unit 104, and thecalculation result stored in the calculation result storage unit 108.The physical property calculation unit 111 outputs an obtainedcalculation result to the calculation result storing unit 112.

The physical property calculation unit 111 also executes the physicalproperty calculation process, when called by the imaginary timeevolution calculation unit 107. Here, the physical property calculationunit 111 executes the physical property calculation process using thecalculation condition data stored in the calculation condition storageunit 102, the initial value stored in the initial value storage unit104, and the calculation result stored in the calculation result storageunit 108. The physical property calculation unit 111 outputs an obtainedcalculation result to the calculation result storing unit 112.

The “physical property calculation process” is a process, of calculatinga physical property such as conductance on the basis of an eigenfunctionand an eigenvalue obtained as a result of numerical calculations. Notethat the following physical properties (1) to (8) may also be calculatedusing the eigenfunction, the eigenvalue, and the like obtained as aresult of calculations.

(1) A charge density spatial distribution may be calculated.

(2) Energy dependence of an electronic density of states may becalculated.

(3) A force acting on an atom may be calculated as an example. Thecalculated force acting on the atom may be utilized to execute amolecular dynamics simulation for tracking atomic motion and execute aprocess of searching for an optimum atomic structure or tracking achemical reaction.

(4) An adiabatic potential curve may be calculated from total energy ofa system, with an elastic modulus being calculated from the calculatedadiabatic potential curve.

(5) An atomic vibrational spectrum may be calculated. Further, aninfrared absorption spectrum may be calculated using the calculatedatomic vibrational spectrum.

(6) An optical physical property such as reflectivity and permittivitymay be calculated.

(7) Excitation energy may be calculated.

(8) A magnetic physical property may be calculated from an electronicspin state.

The calculation result storing unit 112 stores the calculation resultoutputted from the physical property calculation unit 111.

The calculation result output unit 113 outputs the calculation resultstored in the calculation result storage unit 108, the analysis resultstored in the analysis result storage unit 110, and the calculationresult stored in the calculation result storing unit 112, when called bythe real time evolution calculation unit 106. Here, the calculationresult output unit 113 outputs these results to the operator via anoutput device (not illustrated), or outputs these results to a file.

The calculation result output unit 113 also outputs the calculationresult stored in the calculation result storage unit 108 and thecalculation result stored in the calculation result storing unit 112,when called by the imaginary time evolution calculation unit 107. Here,the calculation result output unit 113 outputs these results to theoperator via the output device (not illustrated), or outputs theseresults to a file.

A numerical simulation process executed in the numerical simulationapparatus 100 in this embodiment is described below.

FIG. 3 shows the numerical simulation process executed in the numericalsimulation apparatus 100 in this embodiment. As shown in FIG. 3, thenumerical simulation apparatus 100 sets the calculation condition (StepS101). In detail, the calculation condition setting unit 101 receivesthe calculation condition data from the operator via the input device,or reads the calculation condition data from the file containing thecalculation condition data. The calculation condition setting unit 101outputs the calculation condition data received from the operator or thecalculation condition data read from the file, to the calculationcondition storage unit 102. The calculation condition storage unit 102stores the calculation condition data outputted from the calculationcondition setting unit 101.

Next, the numerical simulation apparatus 100 sets the initial value(Step S102). In detail, the initial value setting unit 103 receives theinitial value from the operator via the input device, or reads theinitial value from the file containing the initial value. The initialvalue setting unit 103 outputs the initial value received from theoperator or the initial value read from the file, to the initial valuestorage unit 104. The initial value storage unit 104 stores the initialvalue outputted from the initial value setting unit 103.

Next, the numerical simulation apparatus 100 judges whether or not toexecute the real time evolution calculation process (Step S103). Indetail, the real time evolution judgment unit 105 receives theinstruction data from the operator via the input device, or reads theinstruction data from the file containing the instruction data. The realtime evolution judgment unit 105 judges whether or not to execute thereal time evolution calculation process, on the basis of the instructiondata received from the operator or the instruction data read from thefile.

Next, when judging to execute the real time evolution calculationprocess (Step S103: Yes), the numerical simulation apparatus 100executes the real time evolution calculation process (Step S104). Indetail, the real time evolution judgment unit 105 calls the real timeevolution calculation unit 106 to execute the real time evolutioncalculation process.

After the real time evolution calculation process ends, the numericalsimulation apparatus 100 further executes the frequency analysis process(Step S105). The numerical simulation apparatus 100 also executes thephysical property calculation process (Step S106). In detail, the realtime evolution calculation unit 106 calls the frequency analysis unit109 to execute the frequency analysis process. The real time evolutioncalculation unit 106 also calls the physical property calculation unit111 to execute the physical property calculation process.

The numerical simulation apparatus 100 then outputs the calculationresults (Step S107). In detail, after the real time evolutioncalculation process, the frequency analysis process, and the physicalproperty calculation process end, the real time evolution calculationunit 106 calls the calculation result output unit 113 to output thecalculation result, the analysis result, and the calculation resultobtained as a result of the series of processes.

When judging not to execute the real time evolution calculation process(Step S103: No), on the other hand, the numerical simulation apparatus100 executes the imaginary time evolution calculation process (StepS108). In detail, the real time evolution judgment unit 105 calls theimaginary time evolution calculation unit 107 to execute the imaginarytime evolution calculation process.

After the imaginary time evolution calculation process ends, thenumerical simulation apparatus 100 further executes the physicalproperty calculation process (Step S106). In detail, the imaginary timeevolution calculation unit 107 calls the physical property calculationunit 111 to execute the physical property calculation process.

The numerical simulation apparatus 100 then outputs the calculationresults (Step S107). In detail, after the imaginary time evolutioncalculation process and the physical property calculation process end,the imaginary time evolution calculation unit 107 calls the calculationresult output unit 113 to output the calculation result and thecalculation result obtained as a result of the series of processes.

Features of the numerical simulation apparatus 100 in this embodimentare described next. The numerical simulation apparatus 100 has thefollowing features (1) to (4).

(1) Numerical calculations of time evolution of a wave functionaccording to the time dependent Schrödinger equation can be executedwith high speed and precision.

(2) Numerical calculations of a steady-state solution of a scatteringproblem can be executed with high speed and precision.

(3) Numerical calculations of imaginary time evolution of a wavefunction used when determining a ground state and an excited state canbe performed with high speed and precision.

(4) A large calculation region needed when setting an infinite boundarycondition can be significantly reduced.

First, the capability of the numerical simulation apparatus 100 in thisembodiment to perform numerical calculations of time evolution of a wavefunction according to the time dependent Schrödinger equation with highspeed and precision is described below.

When performing integration of the above Expression (5), a centraldifference approximation in the real-space finite-difference method isapplied. As a result, as shown in the following Expression (6), the wavefunction after short time period Δt can be approximated using the Besselfunction shown in the following Expression (7). This contributes toreductions in computational complexity of the real time evolutioncalculation process and the imaginary time evolution calculationprocess, with it being possible to execute a numerical simulation withhigh speed and precision.

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 6} \rbrack & \; \\{{\Psi( {r_{j},{t + {\Delta\; t}}} )} = {\sum\limits_{j_{x}^{\prime}}{\sum\limits_{j_{y}^{\prime}}{\sum\limits_{j_{z}^{\prime}}{\lbrack {\prod\limits_{{l = x},y,z}{{\mathbb{i}}^{j_{l} - j_{l}^{\prime}}{\mathbb{e}}^{{- {\mathbb{i}}}\frac{\Delta\; t}{h_{l}^{2}}}{J_{j_{l} - j_{l}^{\prime}}( \frac{\Delta\; t}{h_{l}^{2}} )}}} \rbrack \times {\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({r_{j},t})}}}{\Psi( {r_{j}^{\prime},t} )}}}}}} & (6) \\{\mspace{79mu}\lbrack {{Expression}\mspace{14mu} 7} \rbrack} & \; \\{\mspace{79mu}{{J_{n}(x)} = {\sum\limits_{s = 0}^{\infty}\frac{( {- 1} )^{s}( {x/2} )^{n + {2s}}}{{s!}{( {n + s} )!}}}}} & (7)\end{matrix}$

r_(j) is given by the following Expression (8), and r_(j′) is given bythe following Expression (9). In is a space step size in the l (=x, y,z) direction, and j_(l) and j_(l′) are integers.[Expression 8]r _(j)=(h _(x) j _(x) ,h _(y) j _(y) ,h _(z) j _(z))  (8)[Expression 9]r _(j′)=(h _(x) j′ _(x) ,h _(y) j′ _(y) ,h _(z) j′ _(z))  (9)

FIG. 4 shows a value of the Bessel function J_(j-j′)(Δt/h_(x) ²) and areal part of an integrand exp(i(x_(j)−x_(j′))²/2Δt) of the aboveExpression (5), where x_(j)=jh_(x), x_(j′)=j′h_(x), Δt=0.2, andh_(x)=0.5.

As shown in FIG. 4, the Bessel function J_(j-j′)(x) rapidly decays withincreasing |j−j′| (graph 122). Accordingly, summation can be stopped atsuch a small term that does not influence a calculation result. Thisbeing so, by approximating the above Expression (5) using the aboveExpression (6), a sum relating to j′ of the above Expression (6) iscalculated, and therefore computational complexity can be significantlyreduced when compared with the case of calculating a direct integral. Asa result, the numerical simulation can be executed with high speed.

This is because calculating the direct integral of the above Expression(5) requires integration to be performed over an extremely wide range.Take an x′ component of the integrand of the above Expression (5) givenby the following Expression (10), as one example. In this case, theintegrand for integration of x′ given by the following Expression (11)is a vibrational function (trigonometric function) for x′. This does notdecay with increasing x′, as shown in FIG. 4 (graph 121). Hence, it isnecessary to perform integration over an extremely wide range.[Expression 10]exp[i(r−r′)²/2Δt]  (10)[Expression 11]exp[i(x−x′)²/2Δt]  (11)

Note that an explanation on how the correct result can be obtained usingthe above Expression (6) will be given later.

Next, the capability of the numerical simulation apparatus 100 in thisembodiment to numerically simulate numerical calculations of asteady-state solution of a scattering problem with high speed andprecision is described below. A system in which potential V does notdepend on time is used as an example. In this case, the frequencyanalysis unit 109 can find a steady-state solution of a scatteringproblem using the above Expression (6). Here, a pulse function on planez=z₀ is used as an initial wave function at initial time to, as shown inthe following Expression (12).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 12} \rbrack & \mspace{11mu} \\{{\chi( {r_{j},t_{0}} )} = \{ \begin{matrix}{\frac{1}{\sqrt{N_{x}N_{y}}}{\exp\lbrack {{\mathbb{i}}( {{G_{nx}h_{x}j_{x}} + {G_{ny}h_{y}j_{y}}} )} \rbrack}} & ( {j_{z} = 0} ) \\{0\mspace{365mu}} & ( {j_{z} \neq 0} )\end{matrix} } & (12)\end{matrix}$

G_(nx) and G_(ny) are reciprocal lattice vectors, and are respectivelygiven by the following Expressions (13) and (14). N_(x) is the number oflattice points in the x direction, and n_(x) is an integer. N_(y) is thenumber of lattice points in the y direction, and n_(y) is an integer.

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 13} \rbrack & \; \\{G_{nx} = {\frac{2\pi}{N_{x}h_{x}}n_{x}}} & (13) \\\lbrack {{Expression}\mspace{14mu} 14} \rbrack & \; \\{G_{ny} = {\frac{2\pi}{N_{y}h_{y}}n_{y}}} & (14)\end{matrix}$

The frequency analysis unit 109 can calculate wave function X(r_(j), t)which is evolved into an arbitrary time from initial time t₀, by usingthe initial wave function indicated by the above Expression (12) and theabove Expression (5). Further, the frequency analysis unit 109 canperform a frequency analysis of wave function X(r_(j), t) for energy E,as shown in the following Expression (15).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 15} \rbrack & \; \\{{\Psi( {r_{j};E} )} = {{C(E)}{\int_{t_{0}}^{\infty}{{\mathbb{e}}^{{\mathbb{i}}\;{E{({t - t_{0}})}}}\ {\chi( {r_{j},l} )}{\mathbb{d}l}}}}} & (15)\end{matrix}$

C(E) is given by the following Expression (16). v_(z)(E) is a groupvelocity in the z direction, and is given by the following Expression(17). k_(z) is a wave number in the z direction, and is given by thefollowing Expression (18).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 16} \rbrack & \; \\{{C(E)} = {\frac{1}{h_{z}}{\upsilon_{z}(E)}{\mathbb{e}}^{{\mathbb{i}}\; k_{z}z_{0}}}} & (16) \\\lbrack {{Expression}\mspace{14mu} 17} \rbrack & \; \\{{\upsilon_{z}(E)} = {\frac{1}{h_{z}}\sin\mspace{11mu} k_{z}h_{z}}} & (17) \\\lbrack {{Expression}\mspace{14mu} 18} \rbrack & \; \\{k_{z} = {\frac{1}{h_{z}}{\cos^{- 1}( {1 - {h_{z}^{2}E}} )}}} & (18)\end{matrix}$

The above Expression (15) is used to numerically simulate the wavefunction evolved with time, with the pulse function indicated by theabove Expression (12) being set as the initial wave function. In thisway, steady-state wave function Ψ(r_(j); E) in total energy E can benumerically simulated.

Next, the capability of the numerical simulation apparatus 100 in thisembodiment to numerically simulate numerical calculations of imaginarytime evolution of a wave function used when determining a ground stateand an excited state with high speed and precision is described below. Asystem in which potential V does not depend on time is used as anexample. In this case, the imaginary time evolution calculation unit 107can extend to a calculation method of finding a “steady-state solutionnot dependent on time” in a discrete eigenvalue problem with high speed,by using the above Expression (6) and imaginary time evolutioncalculations. A propagator is defined by the following Expression (19),with eigenfunction Φ_(i)(r) and eigenvalue E_(i) of the system.

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 19} \rbrack & \; \\{{K( {r,{t;r^{\prime}},t_{0}} )} = {\sum\limits_{i}{{\phi_{i}(r)}{\phi_{i}^{*}( r^{\prime} )}{\mathbb{e}}^{{- {\mathbb{i}}}\;{E_{i}{({t - t_{0}})}}}}}} & (19)\end{matrix}$

Integrating for r where r′=r yields the following Expression (20), dueto orthonormality of the eigenfunction.

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 20} \rbrack & \; \\\begin{matrix}{{\int_{- \infty}^{\infty}{{K( {r,{l;r},l_{0}} )}\ {\mathbb{d}^{3}r}}} = {\int_{- \infty}^{\infty}{\sum\limits_{i}{{\phi_{i}(r)}{\phi_{i}^{*}( r^{\prime} )}{\mathbb{e}}^{{- {\mathbb{i}}}\;{E_{i}{({t - t_{0}})}}}\ {\mathbb{d}^{3}r}}}}} \\{= {\sum\limits_{i}{\mathbb{e}}^{{- {\mathbb{i}}}\;{E_{i}{({t - t_{0}})}}}}}\end{matrix} & (20)\end{matrix}$

Furthermore, by formally introducing imaginary time T=it and evolvingover a sufficiently long imaginary time, that is, T→∞, an approximationcan be made as indicated by the following Expression (21). This isbecause the term of lowest energy E₀ decays most slowly.

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 21} \rbrack & \; \\{{\lim\limits_{\tau->\infty}{\int_{- \infty}^{\infty}{{K( {r,{{- {\mathbb{i}\tau}};r},{- {\mathbb{i}\tau}_{0}}} )}\ {\mathbb{d}^{3}r}}}} \approx {\mathbb{e}}^{- {E_{0}{({\tau - \tau_{0}})}}}} & (21)\end{matrix}$

This demonstrates that a discrete eigenstate can be obtained byimaginary time evolution calculations, which is a well-known calculationmethod (for example, see Non-patent Reference 2).

Accordingly, substituting imaginary time for the time in the aboveExpression (6) gives the following Expression (22).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 22} \rbrack & \; \\{{\Psi( {r_{j},{- {{\mathbb{i}}( {\tau + {\Delta\tau}} )}}} )} = {\sum\limits_{j_{x}^{\prime}}{\sum\limits_{j_{y}^{\prime}}{\sum\limits_{j_{z}^{\prime}}{\lbrack {\underset{{l = x},y,z}{\prod}{\mathbb{i}}^{j_{l} - j_{l}^{\prime}}{\mathbb{e}}^{- \frac{\Delta\tau}{h_{l}^{2}}}{J_{j_{l} - j_{l}^{\prime}}( {{- {\mathbb{i}}}\frac{\Delta\tau}{h_{l}^{2}}} )}} \rbrack \times {\mathbb{e}}^{{- {\mathbb{i}\Delta\tau}}\;{V{({r_{j},{- {\mathbb{i}\tau}}})}}}{\Psi( {r_{j}^{\prime};{- {\mathbb{i}\tau}}} )}}}}}} & (22)\end{matrix}$

The imaginary time evolution can be calculated with high speed by usingthe above Expression (22), and also the excited state can be calculatedby using the initial wave function orthogonal to the ground state wavefunction.

Next, the capability of the numerical simulation apparatus 100 in thisembodiment to significantly reduce a large calculation region neededwhen setting an infinite boundary condition is described below.

For example, to obtain a correct result in a scattering problem, asufficiently large space (calculation region) equal to or more thanabout 100 times a spatial size of a scattering potential needs to be setas a calculation condition. This being the case, by applying such an“absorbing boundary condition” that smoothly decays the wave function toan absorbing region, the size of the calculation region can be reduced.

FIG. 5 shows a setting example of a coefficient by which the wavefunction is multiplied in the absorbing boundary condition. As shown inFIG. 5, by applying such an “absorbing boundary condition” that smoothlydecays the wave function by, for example, multiplying the wave functionby the coefficient in absorbing regions 132 a and 132 b, the size of thecalculation region can be reduced to about several times to 10 times thesize of the scatterer. This contributes to reductions in computationalcomplexity of the real time evolution calculation process and theimaginary time evolution calculation process, with it being possible toexecute a numerical simulation with high speed and precision.

Results of executing numerical simulations on simple models using thenumerical simulation apparatus in this embodiment are described below.

FIG. 6 shows a one-dimensional scattering problem by a potentialbarrier. As shown in FIG. 6, suppose an incident wave 142 enters apotential barrier 141 from the left, where time step Δt=0.2 (a.u.), gridwidth h=0.5 (a.u.), size in z direction Lz=1000.0, number of grid pointsin z direction Nz=2000, potential V=0.3 (a.u.), barrier size m=21, andintegral upper limit T=400 (a.u.) (table 145). Moreover, 1 (a.u.)=0.024(fs)=27.211 (eV)=0.52918 Å (text box 146).

FIG. 7 shows a result of analyzing a frequency of the wave functionevolved with time and a result of calculating transmittance to thepotential barrier. As shown in FIG. 7, as a result of analyzing thefrequency of the time-evolved wave function by energy (graphs 148 a, 148b, and 148 c), the transmittance matches the exact solution (graph 147).

FIG. 8 shows a two-dimensional wire model. As shown in FIG. 8, supposean incident wave enters a wire 151 of width Wx and length Wz connectinga left electrode 152 and a right electrode 153, from the left electrode152, where time step Δt=0.2 (a.u.), grid width h=0.5 (a.u.), size in xdirection Lx=4.0, number of grid points in x direction Nx=8, size in zdirection Lz=500.0, number of grid points in z direction Nz=1000,integral upper limit T=200 (a.u.), and incident electron energy is 0.05to 1.10 (table 155). Moreover, 1 (a.u.)=0.024 (fs)=27.211 (eV)=0.52918 Å(text box 156).

FIG. 9 shows a result of analyzing a frequency in the two-dimensionalwire model and a result of calculating conductance. As shown in FIG. 9,by analyzing the frequency of the time-evolved wave function by energy(graphs 158 a to 158 e), the result matching that of the OBM method isobtained (graph 157).

FIG. 10 shows a scattering wave function obtained by applying theabsorbing boundary condition. As shown in FIG. 10, computationalcomplexity can be significantly reduced while maintaining computationalprecision.

FIG. 11 shows calculation cost evaluations. As shown in FIG. 11, giventhat a scattering solution at arbitrary energy can be calculatedaccording to the present method (graph 171), the present method canrealize a faster numerical simulation than the OBM is method (graph172). Which is to say, the present method can numerically simulatehigh-speed electrical conduction properties.

As described above, according to the numerical simulation apparatus 100in this embodiment, by employing the wave function expressed using theBessel function, numerical calculations of time evolution of a wavefunction according to the time dependent Schrödinger equation can beperformed with high speed and precision. In addition, a steady-statesolution of a scattering problem can be obtained with high speed andprecision by numerical calculations. Further, numerical calculations ofimaginary time evolution of a wave function used when determining aground state and an excited state can be performed with high speed andprecision. Moreover, a large calculation region needed when setting aninfinite boundary condition can be significantly reduced.

Furthermore, by applying to the density functional theory which is adominant method of electronic state calculations today, a numericalsolution of the Kohn-Sham equation can be obtained with high speed. Thismakes it possible to numerically simulate dynamic behavior of electronsfor an actual model. Thus, the numerical simulation apparatus 100 can beused for design of nanoscale electronic devices and quantum computers.

A supplementary explanation on the above Expression (6) is given below.

In the path integral representation, the time evolution of the wavefunction after short time period Δt is expressed by the followingExpression (23) (for example, see Non-patent Reference 2).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 23} \rbrack & \; \\{{\Psi( {r,{t + {\Delta\; t}}} )} = {\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; t}}{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{\mathbb{i}}\frac{{({r - r^{\prime}})}^{2}}{2\Delta\; t}}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({\frac{r + r^{\prime}}{2},{t + \frac{\Delta\; t}{2}}})}}}\ {\Psi( {r^{\prime},t} )}{\mathbb{d}^{3}r^{\prime}}}}}} & (23)\end{matrix}$

It is assumed here that potential V depends on time in general. When x,y, and z are used as coordinate variables, the above Expression (23) iswritten as the following Expression (24).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 24} \rbrack & \; \\{{\Psi( {x,y,z,{t + {\Delta\; t}}} )} = {\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; t}}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{\mathbb{i}}\frac{{({x - x^{\prime}})}^{2}}{2\Delta\; t}}{\mathbb{e}}^{{\mathbb{i}}\frac{{({y - y^{\prime}})}^{2}}{2\Delta\; t}}{\mathbb{e}}^{{\mathbb{i}}\frac{{({z - z^{\prime}})}^{2}}{2\Delta\; t}}\  \times {\mathbb{e}}^{{- {\mathbb{i}\Delta}}\; t\;{V{({\frac{x + x^{\prime}}{2},\frac{y + y^{\prime}}{2},\frac{z + z^{\prime}}{2},{t + \frac{\Delta\; t}{2}}})}}}{\Psi( {x^{\prime},y^{\prime},z^{\prime},t} )}{\mathbb{d}x^{\prime}}\ {\mathbb{d}y^{\prime}}\ {\mathbb{d}z^{\prime}}}}}}}} & (24)\end{matrix}$

The following describes only the x coordinate, as this is no loss ofgenerality. That is, the above Expression (24) is simplified to thefollowing Expression (25).

$\begin{matrix}\lbrack {{Expression}\mspace{20mu} 25} \rbrack & \; \\{{\Psi( {x,{t + {\Delta\; t}}} )} = {\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; t}}{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{\mathbb{i}}\frac{{({x - x^{\prime}})}^{2}}{2\Delta\; t}}\ {\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({\frac{x + x^{\prime}}{2},{t + \frac{\Delta\; t}{2}}})}}}{\Psi( {x^{\prime},t} )}{\mathbb{d}x^{\prime}}}}}} & (25)\end{matrix}$

Since contributions to an integral about x′ come only from aneighborhood of x (for example, see Non-patent Reference 2), anintegration variable is transformed from x′ to η such that x′=x+η. As aresult, the above Expression (25) is written as the following Expression(26).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 26} \rbrack & \; \\{{\Psi( {x,{t + {\Delta\; t}}} )} = {\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; t}}{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{lV}{({{x + \frac{\eta}{2}},{l + \frac{\Delta\; t}{2}}})}}}\ {\Psi( {{x + \eta},t} )}{\mathbb{d}\eta}}}}} & (26)\end{matrix}$

Moreover, since main contributions to an integral of η occur only whereη is small (for example, see Non-patent Reference 2), wave functionΨ(x+η, t) is subject to Taylor expansion around x. In addition,potential V(x+η/2, t+Δt/2) is subject to Taylor expansion around t andx. This allows the above Expression (26) to be approximated by thefollowing Expression (27).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 27} \rbrack & \; \\{{\Psi( {x,{t + {\Delta\; t}}} )} \approx {\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; t}}{\int_{- \infty}^{\infty}{{{\mathbb{e}}^{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}}\lbrack {{\Psi( {x,t} )} + {\eta\frac{\partial}{\partial x}{\Psi( {x,t} )}} + {\frac{1}{2}\eta^{2}\frac{\partial^{2}}{\partial x^{2}}{\Psi( {x,t} )}}} \rbrack} \times {\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({{x + \frac{\eta}{2}},t})}}}\ {\mathbb{d}\eta}}}} \approx {\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; t}}{\int_{- \infty}^{\infty}{{{\mathbb{e}}^{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}}\lbrack {{\Psi( {x,t} )} + {\eta\frac{\partial}{\partial x}{\Psi( {x,t} )}} + {\frac{1}{2}\eta^{2}\frac{\partial^{2}}{\partial x^{2}}{\Psi( {x,t} )}} +} \rbrack} \times {\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{t{\lbrack{{V{({x,t})}} + {\frac{\eta}{2}\frac{\partial}{\partial x}{V{({x,t})}}} + {\frac{1}{2}{(\frac{\eta}{2})}^{2}\frac{\partial^{2}}{\partial x^{2}}{V{({x,t})}}}}\rbrack}}}\ {\mathbb{d}\eta}}}}} & (27)\end{matrix}$

Furthermore, a term containing a space derivative of potential V isexpanded using the following Expression (28). This allows the aboveExpression (27) to be approximated by the following Expression (29).[Expression 28]e ^(z)≈1+z  (28)

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 29} \rbrack & \; \\{{\Psi( {x,{t + {\Delta\; t}}} )} \approx {\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; t}}{\int_{- \infty}^{\infty}{{{\mathbb{e}}^{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}}\ \lbrack {{\Psi( {x,t} )} + {\eta\frac{\partial}{\partial x}{\Psi( {x,t} )}} + {\frac{1}{2}\eta^{2}\frac{\partial^{2}}{\partial x^{2}}{\Psi( {x,t} )}}} \rbrack} \times {{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{lV}{({x,t})}}}\lbrack {1 - {{\mathbb{i}\Delta}\;{t( {{\frac{\eta}{2}\frac{\partial}{\partial x}{V( {x,t} )}} + {\frac{1}{2}( \frac{\eta}{2} )^{2}\frac{\partial^{2}}{\partial x^{2}}{V( {x,t} )}}} )}}} \rbrack}{\mathbb{d}\eta}}}}} & (29)\end{matrix}$

Here, components of a product of an even function given by the followingExpression (30) and an odd function given by the following Expression(31) are 0, as shown in the following Expression (32).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 30} \rbrack & \; \\{\mathbb{e}}^{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}} & (30) \\\lbrack {{Expression}\mspace{14mu} 31} \rbrack & \; \\{\eta,\eta^{3},\eta^{5}} & (31) \\\lbrack {{Expression}\mspace{14mu} 32} \rbrack & \; \\{{{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}}\ \eta{\mathbb{d}\eta}}} = 0},{{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}}\ \eta^{3}{\mathbb{d}\eta}}} = 0},{{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}}\ \eta^{5}{\mathbb{d}\eta}}} = 0}} & (32)\end{matrix}$

As shown in the following Expression (33), when up to a first term of Δtis taken into account, and a term containing η⁴ is ignored with η beingregarded as a small quantity, the following Expression (34) is derived.

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 33} \rbrack & \; \\{\eta^{2},{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}}\ \eta^{2}{\mathbb{d}\eta}}},{{\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; t}}{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}}\ \eta^{2}{\mathbb{d}\eta}}}} = {{\mathbb{i}\Delta}\; t}}} & (33) \\\lbrack {{Expression}\mspace{14mu} 34} \rbrack & \; \\{{\Psi( {x,{t + {\Delta\; t}}} )} \approx {\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; t}}{\int_{- \infty}^{\infty}{{{\mathbb{e}}^{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}}\ \lbrack {{\Psi( {x,t} )} + {\frac{1}{2}\eta^{2}\frac{\partial^{2}}{\partial x^{2}}{\Psi( {x,t} )}}} \rbrack}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({x,t})}}}{\mathbb{d}\eta}}}}} & (34)\end{matrix}$

Next, assuming that the space is a discrete space having a periodicboundary condition, the wave function is subject to Fourier seriesexpansion as shown in the following Expression (35), where k_(n) isgiven by the following Expression (36) and x_(j) is given by thefollowing Expression (37).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 35} \rbrack & \; \\{{\Psi( {x_{j},t} )} = {\sum\limits_{n = {- \frac{N_{x}}{2}}}^{\frac{N_{x}}{2} - 1}{c_{n}{\mathbb{e}}^{{{\mathbb{i}}k}_{n}x_{j}}}}} & (35) \\\lbrack {{Expression}\mspace{14mu} 36} \rbrack & \; \\{k_{n} = \frac{2n\;\pi}{L_{x}}} & (36) \\\lbrack {{Expression}\mspace{14mu} 37} \rbrack & \; \\{x_{j} = {h_{x}{j( {{j = {- \frac{N_{x}}{2}}},\ldots\mspace{14mu},{- 1},0,1,\ldots\mspace{14mu},{\frac{N_{x}}{2} - 1}} )}}} & (37)\end{matrix}$

c_(n) is an expansion coefficient. L_(x) is the size in the x direction,where L_(x)=N_(x)h_(x). Subsequently, the size of the space grows atinfinity, such that L_(x)→∞. For simplicity, N_(x) is assumed to be aneven number.

Further, applying the central difference approximation in the real-spacefinite-difference method to ∂²Ψ(x, t)/∂x² yields the followingExpression (38). From this, the time evolution of the wave function inthe discrete space is indicated by the following Expression (39).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 38} \rbrack & \; \\\begin{matrix}{{\frac{\partial^{2}}{\partial x^{2}}{\Psi( {x_{j},t} )}} \approx \frac{{\Psi( {{x_{j} - h_{x}},t} )} - {2{\Psi( {x_{j},t} )}} + {\Psi( {{x_{j} + h_{x}},t} )}}{h_{x}^{2}}} \\{= {\sum\limits_{n}{c_{n}{\frac{1}{h_{x}^{2}}\lbrack {{\mathbb{e}}^{{{\mathbb{i}}\;{k_{n}{({x_{j} - h_{x}})}}} -}2{\mathbb{e}}^{{{\mathbb{i}}\; k_{n}x_{j}} +}{\mathbb{e}}^{{\mathbb{i}}\;{k_{n}{({x_{j} + h_{x}})}}}} \rbrack}}}} \\{= {\sum\limits_{n}{c_{n}\frac{1}{h_{x}^{2}}{{\mathbb{e}}^{{\mathbb{i}}\; k_{n}x_{j}}\lbrack {{\mathbb{e}}^{{\mathbb{i}}\; k_{n}h_{x}} - 2 + {\mathbb{e}}^{{\mathbb{i}}\;{k\;}_{n}h_{x}}} \rbrack}}}} \\{= {\sum\limits_{n}{c_{n}\frac{- 2}{h_{x}^{2}}{{\mathbb{e}}^{{\mathbb{i}}\; k_{n}x_{j}}\lbrack {1 - {\cos( {k_{n}h_{x}} )}} \rbrack}}}}\end{matrix} & (38) \\\lbrack {{Expression}\mspace{14mu} 39} \rbrack & \; \\{{\Psi( {x_{j},{t + {\Delta\; t}}} )} \approx {\sum\limits_{n}{\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; l}}{\int_{- \infty}^{\infty}{{{\mathbb{e}}^{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}}\ \lbrack {1 - {\frac{\eta^{2}}{h_{x}^{2}}\lbrack {1 - {\cos( {k_{n}h_{x}} )}} \rbrack}} \rbrack} \times c_{n}{\mathbb{e}}^{{\mathbb{i}}\; k_{n}x_{j}}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({x_{j},t})}}}{\mathbb{d}\eta}}}}}} & (39)\end{matrix}$

The use of the above Expression (28) allows the above Expression (39) tobe approximated by the following Expression (40).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 40} \rbrack & \; \\\begin{matrix}{{{\Psi( {x_{j},{t + {\Delta\; t}}} )} \approx {\sum\limits_{n}{\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; t}}{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}}{\mathbb{e}}^{- {\frac{\eta^{2}}{h_{x}^{2}}{\lbrack{1 - {\cos{({k_{n}h_{x}})}}}\rbrack}}}}}}}}\ } \\{c_{n}{\mathbb{e}}^{{{\mathbb{i}}k}_{n}x_{j}}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({x_{j},t})}}}{\mathbb{d}\eta}} \\{{= {\sum\limits_{n}{\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; t}}{\int_{- \infty}^{\infty}{\mathbb{e}}^{{{\mathbb{i}}\frac{\eta^{2}}{2\Delta\; t}} - {\frac{\eta^{2}}{h_{x}^{2}}{\lbrack{1 - {\cos{({h_{n}h_{x}})}}}\rbrack}}}}}}}\ } \\{c_{n}{\mathbb{e}}^{{{\mathbb{i}}k}_{n}x_{j}}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\; l\;{V{({x_{j},l})}}}{\mathbb{d}\eta}} \\{= {\sum\limits_{n}{\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; t}}{\int_{- \infty}^{\infty}{\mathbb{e}}^{- {{\mathbb{i}\eta}^{2}{\lbrack{- {\frac{1}{2\Delta\; t}{\lbrack{1 + {2{\mathbb{i}\Delta}\; t\;{\frac{1}{h_{x}^{2}}{\lbrack{1 - {\cos{({k_{n}h_{x}})}}}\rbrack}}}}\rbrack}}}\rbrack}}}}}}} \\{c_{n}{\mathbb{e}}^{{\mathbb{i}}\; k_{n}x_{j}}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{lV}{({x_{j},l})}}}{\mathbb{d}\eta}}\end{matrix} & (40)\end{matrix}$

By performing integration about η according to an integral formuladefined by the following Expression (41), the above Expression (40) canbe approximated by the following Expression (42).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 41} \rbrack & \; \\{{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{- {\mathbb{i}}}\; a\;\eta^{2}}{\mathbb{d}\eta}}} = \sqrt{\frac{\pi}{{\mathbb{i}}\; a}}} & (41) \\\lbrack {{Expression}\mspace{14mu} 42} \rbrack & \; \\{{\Psi( {x_{j},{t + {\Delta\; t}}} )} \approx {\sum\limits_{n}{\sqrt{\frac{1}{2{\pi\mathbb{i}\Delta}\; t}}{\sqrt{\frac{{- 2}{\pi\Delta}\; t}{\mathbb{i}}}\lbrack {1 + {2{\mathbb{i}\Delta}\; t\;{\frac{1}{h_{x}^{2}}\lbrack {1 - {\cos( {k_{n}h_{x}} )}} \rbrack}}} \rbrack}^{- \frac{1}{2}} \times c_{n}{\mathbb{e}}^{{\mathbb{i}}\; k_{n}x_{j}}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({x_{j},t})}}}}}} & (42)\end{matrix}$

Moreover, since Δt is small, by using the above Expression (28) again,the above Expression (42) can be approximated by the followingExpression (43).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 43} \rbrack & \; \\{{\Psi( {x_{j},{t + {\Delta\; t}}} )} \approx {\sum\limits_{n}{{\mathbb{e}}^{{- {\mathbb{i}}}\;\Delta\; t\;{\frac{1}{h_{x}^{2}}{\lbrack{1 - {\cos{({k_{n}h_{x}})}}}\rbrack}}}c_{n}{\mathbb{e}}^{{\mathbb{i}}\; k_{n}x_{j}}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({x_{j},l})}}}}}} & (43)\end{matrix}$

Fourier expansion coefficient c_(n) is given by the following Expression(44). Accordingly, the above Expression (43) is indicated by thefollowing Expression (45).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 44} \rbrack & \; \\{c_{n} = {\frac{1}{N_{x}}{\sum\limits_{j^{\prime} = \frac{N_{x}}{2}}^{\frac{N_{x}}{2} - 1}{{\Psi( {x_{j},t} )}{\mathbb{e}}^{{- {\mathbb{i}}}\;\frac{2n\;\pi}{L_{x}}x_{j^{\prime}}}}}}} & (44) \\\lbrack {{Expression}\mspace{14mu} 45} \rbrack & \; \\\begin{matrix}{{\Psi( {x_{j},{t + {\Delta\; t}}} )} = {\frac{1}{N_{x}}{\sum\limits_{n}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\; t{\frac{1}{h_{x}^{2}}{\lbrack{1 - {\cos{({k_{n}h_{x}})}}}\rbrack}}}}}} \\{\sum\limits_{j^{\prime}}{{\Psi( {x_{j^{\prime}},t} )}{\mathbb{e}}^{{- {\mathbb{i}}}\frac{2\; n\;\pi}{L_{x}}{x_{j}}^{\prime}}{\mathbb{e}}^{{\mathbb{i}}\frac{2n\;\pi}{L_{x}}x_{j}}\;{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({x_{j},t})}}}}} \\{= {\frac{1}{N_{x}}{\sum\limits_{j^{\prime}}{\sum\limits_{n}{{\mathbb{e}}^{i\frac{2n\;\pi}{L_{x}}{({x_{j} - x_{j^{\prime}}})}}{\mathbb{e}}^{{- i}\;\Delta\; t{\frac{1}{h_{x}^{2}}{\lbrack{1 - {\cos{({k_{n}h_{x}})}}}\rbrack}}}}}}}} \\{{\mathbb{e}}^{{- i}\;\Delta\;{{tV}{({x_{j,}t})}}}{\Psi( {x_{j^{\prime}},t} )}}\end{matrix} & (45)\end{matrix}$

Furthermore, by fixing h_(x) and setting N_(x)→∞, that is,L_(x)=h_(x)N_(x)→∞, to change to representation in a discrete infinitespace, the following Expression (46) is obtained. Here, the followingExpressions (47) to (50) are taken into account.

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 46} \rbrack & \; \\{{\Psi( {x_{j},{t + {\Delta\; t}}} )} = {\sum\limits_{j^{\prime}}{\frac{h_{x}}{2\pi}{\int_{- \frac{\pi}{h_{x}}}^{\frac{\pi}{h_{x}}}{{\mathbb{e}}^{{\mathbb{i}}\;{k{({x_{j} - x_{j^{\prime}}})}}}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\; t{\frac{1}{h_{x}^{2}}{\lbrack{1 - {\cos{({k_{n}h_{x}})}}}\rbrack}}}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({x_{j},t})}}}{\Psi( {x_{j^{\prime}},t} )}{\mathbb{d}k}}}}}} & (46) \\\lbrack {{Expression}\mspace{14mu} 47} \rbrack & \; \\{\mspace{79mu}{\frac{2\pi\; n}{N_{x}h_{x}}->k}} & (47) \\\lbrack {{Expression}\mspace{14mu} 48} \rbrack & \; \\{\mspace{79mu}{\frac{2\pi}{N_{x}h_{x}}->{\mathbb{d}k}}} & (48) \\\lbrack {{Expression}\mspace{14mu} 49} \rbrack & \; \\{\mspace{79mu}{\sum\limits_{n = {- \frac{N_{x}}{2}}}^{\frac{N_{x}}{2} - 1}{->\int_{- \frac{\pi}{h_{x}}}^{\frac{\pi}{h_{x}}}}}} & (49) \\\lbrack {{Expression}\mspace{14mu} 50} \rbrack & \; \\\begin{matrix}{\mspace{79mu}{\frac{1}{N_{x}} = {\frac{h_{x}}{2\pi}\frac{2\pi}{N_{x}h_{x}x}}}} \\{= {\frac{h_{x}}{2\pi}{\mathbb{d}k}}}\end{matrix} & (50)\end{matrix}$

As a result of further performing integration variable transformationθ=h_(x)k, the above Expression (46) is indicated by the followingExpression (51).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 51} \rbrack & \; \\\begin{matrix}{{\Psi( {x_{j},{t + {\Delta\; t}}} )} = {\sum\limits_{j^{\prime}}{\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{\mathbb{e}}^{{\mathbb{i}\theta}{({j - j^{\prime}})}}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\; t\;\frac{1}{h_{x}^{2}}{({1 - {\cos\;\theta}})}}}}}}} \\{{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({x_{j},t})}}}{\Psi( {x_{j^{\prime}},t} )}{\mathbb{d}\theta}} \\{= {\sum\limits_{j^{\prime}}{\frac{1}{\pi}{\int_{0}^{\pi}{{\cos( {j - j^{\prime}} )}{\theta\mathbb{e}}^{{- {\mathbb{i}\Delta}}\; t\;\frac{1}{h_{x}^{2}}{({1 - {\cos\;\theta}})}}}}}}} \\{{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({x_{j},t})}}}{\Psi( {x_{j^{\prime}},t} )}{\mathbb{d}\theta}}\end{matrix} & (51)\end{matrix}$

By employing Hansen's integral representation of the Bessel function asshown in the following Expression (52), the above Expression (51) isindicated by the following Expression (53).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 52} \rbrack & \; \\{{J_{n}(z)} = {\frac{1}{{\pi\mathbb{i}}^{n}}{\int_{0}^{\pi}{\lbrack {\cos\; n\;\theta} \rbrack{\mathbb{e}}^{{\mathbb{i}}\; z\;\cos\;\theta}}}}} & (52) \\\lbrack {{Expression}\mspace{14mu} 53} \rbrack & \; \\{{\Psi( {x_{j},{t + {\Delta\; t}}} )} = {\sum\limits_{j^{\prime}}{{\mathbb{i}}^{j - j^{\prime}}{\mathbb{e}}^{{- {\mathbb{i}}}\frac{\Delta\; t}{h_{x}^{2}}}{J_{j - j^{\prime}}( \frac{\Delta\; t}{h_{x}^{2}} )}{\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({x_{j},t})}}}{\Psi( {x_{j^{\prime}},t} )}}}} & (53)\end{matrix}$

J_(n) is the Bessel function, and is given by the above Expression (7).Integration about y and z is the same as above. In the three-dimensionalspace, the above Expression (23) is represented by the followingExpression (54).

$\begin{matrix}\lbrack {{Expression}\mspace{14mu} 54} \rbrack & \; \\{{\Psi( {r_{j},{t + {\Delta\; t}}} )} = {\sum\limits_{j_{x}^{\prime}}{\sum\limits_{j_{y}^{\prime}}{\sum\limits_{j_{z}^{\prime}}{\lbrack {\prod\limits_{{l = x},y,z}{{\mathbb{i}}^{j_{l} - j_{l}^{\prime}}{\mathbb{e}}^{{- {\mathbb{i}}}\frac{\Delta\; t}{h_{l}^{2}}}{J_{j_{l} - j_{l}^{\prime}}( \frac{\Delta\; t}{h_{l}^{2}} )}}} \rbrack \times {\mathbb{e}}^{{- {\mathbb{i}\Delta}}\;{{tV}{({r_{j},t})}}}{\Psi( {r_{j}^{\prime},t} )}}}}}} & (54)\end{matrix}$

Note that the present invention can also apply the following differenceapproximation in the real-space finite-difference method when executingthe integration of the above Expression (5).

$\begin{matrix}{\lbrack {{Expression}\mspace{14mu} 55} \rbrack{{\frac{\partial^{2}}{\partial x^{2}}{\psi( {x_{j},t} )}} = {{\sum\limits_{n = {- N_{f}}}^{N_{f}}{c_{n}{\psi( {{x_{j} + {nh}_{x}},t} )}}} + {O( ( h_{x} )^{{2N} + 2} )}}}} & \;\end{matrix}$

N_(f) is an approximation order, and c_(n) is a weighting coefficient.

In particular, in the case of a simplest central differenceapproximation, the following Expression (56) applies.

[Expression  56]${\frac{\partial^{2}}{\partial x^{2}}{\psi( {x_{j},t} )}} \approx \frac{{\psi( {x_{j + 1},t} )} - {2{\psi( {x_{j},t} )}} + {\psi( {x_{j - 1},t} )}}{h_{x}^{2}}$

The wave function after short time period Δt can be written as thefollowing Expression (57) using the Bessel function (for example, seeNon-patent Reference 3 (Chelikowsky J R, Troullier N, Wu K and Saad Y,Phys. Rev. B, 50(16), 11355 (1994)).

[Expression  57]${\psi( {r_{j},{t + {\Delta\; t}}} )} = {{\mathbb{e}}^{{- {\mathbb{i}}}\;\frac{\Delta\; t}{2}{V{({r_{j},{t + \frac{\Delta\; t}{2}}})}}}{\sum\limits_{j_{x}^{\prime}}{\sum\limits_{j_{y}^{\prime}}{\sum\limits_{j_{z}^{\prime}}{\lbrack {\prod\limits_{{v = x},y,z}{{\mathbb{i}}^{j_{v} - j_{v}^{\prime}}{\mathbb{e}}^{{- {\mathbb{i}}}\;\frac{\Delta\; t}{h_{v}^{2}}}{J_{j_{v} - j_{v}^{\prime}}( \frac{\Delta\; t}{h_{v}^{2}} )}}} \rbrack \times {\mathbb{e}}^{{- {\mathbb{i}}}\;\frac{\Delta\; t}{2}{V{({r_{j^{\prime}},{t + \frac{\Delta\; t}{2}}})}}}{\psi( {r_{j^{\prime}},t} )}}}}}}$

Expression (58) in the above Expression (57) is the same Bessel functionas the above Expression (7). In addition, the same condition ofExpression (59) as the above Expressions (8) and (9) is satisfied.

[Expression  58]${J_{n}(x)} = {\sum\limits_{s = 0}^{\infty}{\frac{( {- 1} )^{s}( {x/2} )^{n + {2s}}}{{s!}{( {n + s} )!}}\lbrack {{Expression}\mspace{14mu} 59} \rbrack}}$r_(j) = (h_(x)j_(x), h_(y)j_(y), h_(z)j_(z)), r_(j^(′)) = (h_(x)j_(x)^(′), h_(y)j_(y)^(′), h_(z)j_(z)^(′))

where h_(v) is a space step size in the v direction, and j_(v) andj_(v)′ are integers.

Meanwhile, in the case of using a difference approximation of the orderN_(f), the following Expression (60) applies.

[Expression  60]${\psi( {r_{j},{t + {\Delta\; t}}} )} = {\frac{1}{\pi^{3}}{\mathbb{e}}^{{- {\mathbb{i}}}\;\frac{\Delta\; t}{2}{V{({r_{j},{t + \frac{\Delta\; t}{2}}})}}} \times {\sum\limits_{j_{x}^{\prime}}{\sum\limits_{j_{y}^{\prime}}{\sum\limits_{j_{z}^{\prime}}{\lbrack {\underset{{v = x},y,z}{\prod}{\int_{0}^{\pi}{\lbrack {{\cos( {j_{v} - j_{v}^{\prime}} )}\theta} \rbrack{\exp( {{\mathbb{i}}\;\frac{\Delta\; t}{h_{v}^{2}}{\sum\limits_{l = 0}^{N_{f}}{{c_{l}( {1 - {\frac{1}{2}\delta_{l0}}} )}{\cos( {l\;\theta} )}}}} )}{\mathbb{d}\theta}}}} \rbrack \times {\quad{{\mathbb{e}}^{{- {\mathbb{i}}}\;\frac{\Delta\; t}{2}{V{({r_{j^{\prime}},{t + \frac{\Delta\; t}{2}}})}}}{\psi( {r_{j^{\prime}},t} )}}}}}}}}$

As described above, Bessel function J_(j-j′)(x) rapidly decays withincreasing |j−j′| as shown in FIG. 4, so that summation can be stoppedat such a small term that does not influence a calculation result. Thiscontributes to a significant reduction in computational complexity forthe integration of Expression (5), that is, for obtaining a sum about J′in Expression (57). Thus, numerical calculations can be performed withhigh speed.

On the other hand, consider the case of the direct integration ofExpression (5). The integrand of Expression (5) is given by Expression(61).[Expression 61]exp[i(r−r′)²/2Δt]

This being so, for example the integrand for the integration of x′ isgiven by Expression (62).[Expression 62]exp[i(x−x′)²/2Δt]

The above Expression (62) is a vibrational function (trigonometricfunction) for x′. As can be understood from its real part shown in FIG.4, there is no decay with increasing x′, and so integration needs to beperformed over an extremely wide range. Therefore, the correct resultcan be obtained by using Expression (57), in the same way as describedabove.

In addition, in the present invention it is possible to extend to acalculation method of finding a “steady-state solution not dependent ontime” in a discrete eigenvalue problem with high speed, by applying theconcept of “imaginary time evolution” to Expression (57). A timeevolution propagator can be written as Expression (63) same as the aboveExpression (19), using eigenfunction Φ_(i)(r) and eigenvalue E_(i) of asystem in question.

$\begin{matrix}{\lbrack {{Expression}\mspace{14mu} 63} \rbrack{{K( {r,{t;r^{\prime}},t_{0}} )} = {\sum\limits_{i}{{\phi_{i}(r)}{\phi_{i}^{*}( r^{\prime} )}{\mathbb{e}}^{{- {\mathbb{i}}}\;{E_{i}{({t - t_{0}})}}}}}}} & \;\end{matrix}$

Integrating for r where r′=r yields Expression (64) same as the aboveExpression (20), due to orthonormality of the eigenfunction.

[Expression  64] $\begin{matrix}{{\int_{- \infty}^{\infty}{{\mathbb{d}^{3}r}\; K( {r,{t;r},t_{0}} )}} = {\int_{- \infty}^{\infty}{{\mathbb{d}^{3}r}{\sum\limits_{i}{{\phi_{i}(r)}{\phi_{i}^{*}(r)}{\mathbb{e}}^{{- {\mathbb{i}}}\;{E_{i}{({t - t_{0}})}}}}}}}} \\{= {\sum\limits_{i}{\mathbb{e}}^{{- {\mathbb{i}}}\;{E_{i}{({t - t_{0}})}}}}}\end{matrix}$

Furthermore, by formally introducing imaginary time T=it and evolvingover a sufficiently long imaginary time, that is, T→∞, Expression (65)same as the above Expression (21) is obtained because the term of lowestenergy E₀ decays most slowly.

[Expression  65]${\lim\limits_{\tau->\infty}{\int_{- \infty}^{\infty}{\mathbb{d}^{3}{{rK}( {r,{{- {\mathbb{i}\tau}};r},{- {\mathbb{i}\tau}_{0}}} )}}}} \approx {\mathbb{e}}^{- {E_{0}{({\tau - \tau_{0}})}}}$

This demonstrates that a discrete eigenstate can be obtained byimaginary time evolution calculations, which is a well-known calculationmethod.

Applying the numerical calculation method of Expression (57) to theabove concept leads to Expression (66).

[Expression  66]${\psi( {r_{j},{- {{\mathbb{i}}( {\tau + {\Delta\tau}} )}}} )} = {{\mathbb{e}}^{{- \frac{\Delta\;\tau}{2}}{V{({r_{j},{t + \frac{\Delta\; t}{2}}})}}}{\sum\limits_{j_{x}^{\prime}}{\sum\limits_{j_{y}^{\prime}}{\sum\limits_{j_{z}^{\prime}}{\lbrack {\prod\limits_{{v = x},y,z}{{\mathbb{i}}^{j_{v} - j_{v}^{\prime}}{\mathbb{e}}^{- \frac{\Delta\; t}{h_{v}^{2}}}{J_{j_{v} - j_{v}^{\prime}}( {{- {\mathbb{i}}}\;\frac{\Delta\tau}{h_{v}^{2}}} )}}} \rbrack \times {\mathbb{e}}^{{- \frac{\Delta\;\tau}{2}}{V{({r_{j^{\prime}},{t + \frac{\Delta\; t}{2}}})}}}{\psi( {r_{j^{\prime}},{- {\mathbb{i}\tau}}} )}}}}}}$

Furthermore, by using Expression (67) about a relation between Besselfunction J_(n)(x) and modified Bessel function I_(n)(x), Expression (68)is obtained.

[Expression  67]      J_(n)(−𝕚 x) = (−𝕚)^(n)I_(n)(x)[Expression  68]${\psi( {r_{j},{- {{\mathbb{i}}( {\tau + {\Delta\tau}} )}}} )} = {{\mathbb{e}}^{{- \frac{\Delta\;\tau}{2}}{V{({r_{j},{t + \frac{\Delta\; t}{2}}})}}} \times {\sum\limits_{j_{x}^{\prime}}{\sum\limits_{j_{y}^{\prime}}{\sum\limits_{j_{z}^{\prime}}{\lbrack {\prod\limits_{{v = x},y,z}{{\mathbb{e}}^{- \frac{\Delta\;\tau}{h_{v}^{2}}}{I_{j_{v} - j_{v}^{\prime}}( \frac{\Delta\tau}{h_{v}^{2}} )}}} \rbrack{\mathbb{e}}^{{- \frac{\Delta\;\tau}{2}}{V{({r_{j^{\prime}},{t + \frac{\Delta\; t}{2}}})}}}{\psi( {r_{j^{\prime}},{- {\mathbb{i}\tau}}} )}}}}}}$

This makes it possible to calculate imaginary time evolution with highspeed, so that the ground state wave function can be obtained with highspeed and precision. Moreover, the excited state can be calculated aswell, by using the initial wave function orthogonal to the ground statewave function.

On the other hand, in the case of using the high-order differenceapproximation shown in the above Expressions (55) to (62), the followingExpression (69) applies.

[Expression  69]${\psi( {r_{j},{- {{\mathbb{i}}( {\tau + {\Delta\tau}} )}}} )} = {\quad{\frac{1}{\pi^{3}}{\mathbb{e}}^{{- \frac{\Delta\;\tau}{2}}{V{({r_{j},{t + \frac{\Delta\; t}{2}}})}}} \times {\sum\limits_{j_{x}^{\prime}}{\sum\limits_{j_{y}^{\prime}}{\sum\limits_{j_{z}^{\prime}}{\lbrack {\prod\limits_{{v = x},y,z}{\int_{0}^{\pi}{\lbrack {{\cos( {j_{v} - j_{v}^{\prime}} )}\theta} \rbrack{\exp( {\frac{\Delta\tau}{h_{v}^{2}}{\sum\limits_{l = 0}^{N_{f}}{( {1 - {\frac{1}{2}\delta_{l0}}} )c_{l}{\cos( {l\;\theta} )}}}} )}{\mathbb{d}\theta}}}} \rbrack \times {\mathbb{e}}^{{- \frac{\Delta\tau}{2}}{V{({r_{j^{\prime}},{t + \frac{\Delta\; t}{2}}})}}}{\psi( {r_{j^{\prime}},{- {\mathbb{i}\tau}}} )}}}}}}}$

Other Variations

Note that the numerical simulation apparatus 100 may be provided with aCentral Processing Unit (CPU), a Random Access Memory (RAM), a Read OnlyMemory (ROM), a Hard Disk Drive (HDD), a network adaptor, and the like.A program for controlling the numerical simulation apparatus 100(hereafter referred to as a numerical simulation program) may beinstalled on the HDD or the like, so that each function of the numericalsimulation apparatus 100 can be realized by executing the numericalsimulation program.

Moreover, the numerical simulation apparatus 100 may be constituted by aplurality of computer systems. In such a case, the plurality of computersystems may be used as one composition computer system, like gridcomputing.

Moreover, the numerical simulation apparatus 100 may be composed of afront end and a back end. In such a case, the initial value, thecalculation condition, and the instruction data are received from theoperator in the front end, and the real time evolution calculationprocess, the imaginary time evolution calculation process, the frequencyanalysis process, the physical property calculation process, and thelike are executed using the received initial value, calculationcondition, and instruction data in the back end.

Moreover, the numerical simulation program may be recorded on arecording medium that is readable by a hardware system such as acomputer system or an embedded system. The numerical simulation programmay also be read and executed on another hardware system via therecording medium, as a result of which each function of the numericalsimulation apparatus 100 can be realized on another hardware system.Examples of the recording medium readable on a computer system includean optical recording medium (for example, a CD-ROM), a magneticrecording medium (for example, a hard disk), a magneto-optical recordingmedium (for example, an MO), and a semiconductor memory (for example, amemory card).

Moreover, the numerical simulation program may be held in a hardwaresystem connected to a network such as an internet, a local area network,or the like. The numerical simulation program may also be downloaded andexecuted on another hardware system via the network, as a result ofwhich each function of the numerical simulation apparatus 100 can berealized on another hardware system. Examples of the network include aterrestrial broadcasting network, a satellite broadcasting network,Power Line Communication (PLC), a mobile telephone network, a wiredcommunication network (for example, IEEE 802.3), and a wirelesscommunication network (for example, IEEE 802.11).

Moreover, the numerical simulation apparatus 100 may be realized by anLSI that includes each function of the numerical simulation apparatus100.

The LSI may be realized by a full-custom Large Scale Integration (LSI),a semi-custom LSI such as an Application Specific Integrated Circuit(ASIC), a programmable logic device such as a FPGA or a CPLD, or adynamic reconfigurable device that enables a circuit structure to bechanged dynamically.

Moreover, design data for implementing each processor function on theLSI may be a program (hereafter referred to as an HDL program) writtenin a hardware description language. The design data may also be agate-level netlist obtained by logical synthesis of the HDL program. Thedesign data may also be macro cell information in which placementinformation, process conditions, and the like are added to thegate-level netlist. The design data may also be mask data specifyingdimensions, timings, and the like. Examples of the hardware descriptionlanguage include Very high speed integrated circuit Hardware DescriptionLanguage (VHDL), Verilog-HDL, SystemC, and the like.

Moreover, the design data may be recorded on a recording medium that isreadable by a hardware system such as a computer system or an embeddedsystem. The design data may also be read and executed on anotherhardware system via the recording medium. The design data read onanother hardware system via the recording medium may be downloaded to aprogrammable logic device via a download cable. Examples of therecording medium readable on a computer system include an opticalrecording medium (for example, a CD-ROM), a magnetic recording medium(for example, a hard disk), a magneto-optical recording medium (forexample, an MO), and a semiconductor memory (for example, a memorycard).

The design data may also be held in a hardware system connected to anetwork such as an internet, a local area network, or the like. Thedesign data may also be downloaded and executed on another hardwaresystem via the network. The design data acquired on another hardwaresystem via the network may be downloaded to a programmable logic devicevia a download cable. Examples of the network include a terrestrialbroadcasting network, a satellite broadcasting network, PLC, a mobiletelephone network, a wired communication network (for example, IEEE802.3), and a wireless communication network (for example, IEEE 802.11).

The design data may also be recorded on a serial ROM so as to betransferred to an FPGA upon activation. The design data recorded on theserial ROM may be directly downloaded to the FPGA upon activation.

Alternatively, the design data may be generated by a microprocessor anddownloaded to the FPGA upon activation.

INDUSTRIAL APPLICABILITY

The present invention can be used as a numerical simulation apparatusand the like needed for developing electronic devices or quantumcomputers in consideration of nanoscale properties, and in particular asa numerical simulation apparatus and the like capable of executing anumerical simulation with high speed and precision by reducingcomputational complexity.

1. A numerical simulation apparatus that executes a numerical simulationusing a second wave function which is a solution of a time dependentSchrödinger equation, said numerical simulation apparatus comprising: atime evolution calculation unit configured to calculate a result of thesecond wave function at a plurality of times from an initial time inincrements of a predetermined time period according to a time evolutioncalculation of the second wave function, wherein the result of thesecond wave function at each time is obtained by applying a centraldifference approximation in a real-space finite-difference method to afirst wave function expressed using a propagator, and the time evolutioncalculation is performed using a Bessel function; and a calculationresult storage unit configured to store the result of the second wavefunction obtained at each time by said time evolution calculation unitwhile evolving the second wave function in increments of thepredetermined time period.
 2. The numerical simulation apparatusaccording to claim 1, further comprising: an initial wave functionstorage unit configured to store a pulse function as the second wavefunction at the initial time; and a frequency analysis unit configuredto calculate a steady-state solution of a scattering problem forpredetermined energy, by performing a frequency analysis of the secondwave function for the predetermined energy using the calculation resultof the second wave function obtained by said time evolution calculationunit from the pulse function stored in said initial wave functionstorage unit.
 3. The numerical simulation apparatus according to claim1, wherein said time evolution calculation unit is configured tocalculate imaginary time evolution of the second wave function, usingimaginary time that is obtained by multiplying time by an imaginarynumber.
 4. The numerical simulation apparatus according to claim 1,further comprising an absorbing boundary condition storage unitconfigured to store an absorption coefficient weighted so as to absorban influence that is exerted on a model by the wave function in aboundary region of the model, the model being subject to the numericalsimulation, wherein said time evolution calculation unit is configuredto calculate the time evolution of the second wave function using theabsorption coefficient stored in said absorbing boundary conditionstorage unit.
 5. A numerical simulation method for executing a numericalsimulation using a second wave function which is a solution of a timedependent Schrödinger equation, said numerical simulation methodcomprising: calculating a result of the second wave function at aplurality of times from an initial time in increments of a predeterminedtime period according to a time evolution calculation of the second wavefunction, wherein the result of the second wave function at each time isobtained by applying a central difference approximation in a real-spacefinite-difference method to a first wave function expressed using apropagator, and the time evolution calculation is performed using aBessel function; and storing the result of the second wave functionobtained at each time in said calculating while evolving the second wavefunction in increments of the predetermined time period.
 6. Anon-transitory computer readable recording medium having stored thereona numerical simulation program for causing a computer system to executea numerical simulation using a second wave function which is a solutionof a time dependent Schrödinger equation, wherein, when executed, saidnumerical simulation program causes the computer system to execute amethod comprising: calculating a result of the second wave function at aplurality of times from an initial time in increments of a predeterminedtime period according to a time evolution calculation of the second wavefunction, wherein the result of the second wave function at each time isobtained by applying a central difference approximation in a real-spacefinite-difference method to a first wave function expressed using apropagator, and the time evolution calculation is performed using aBessel function; and storing the result of the second wave functionobtained at each time in said calculating while evolving the second wavefunction in increments of the predetermined time period.